% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function h = base_plot_q(base,scale)
% h = base_plot_q(base,scale)
%
% Plot the dual basis of base

 meshg = load('plot_mesh.mat');
 % increase the resolution
 pgn = zeros(2,size(meshg.tg,2)*size(meshg.pg,2));
 for i = 1:size(meshg.tg,2)
   x = meshg.pg(1,meshg.tg(1:3,i));
   y = meshg.pg(2,meshg.tg(1:3,i));
   BK = [x(2)-x(1) , x(3)-x(1) ;
         y(2)-y(1) , y(3)-y(1) ];
   XY = BK*meshg.pg;
   XY(1,:) = XY(1,:) + x(1);
   XY(2,:) = XY(2,:) + y(1);
   pgn(:,(i-1)*size(meshg.pg,2)+1 : size(meshg.pg,2)*i) = XY;
 end

 OMEGA = zeros(2,base.mk,size(meshg.pg,2));
 for i = 1:base.mk
   ox_si = base.o_s{1,i};
   oy_si = base.o_s{2,i};
   for j = 1:size(ox_si,2)
     OMEGA(1,i,:) = OMEGA(1,i,:) + ...
                    permute( ...
                    ox_si(1,j)*(meshg.pg(1,:).^ox_si(2,j)) .* ...
                               (meshg.pg(2,:).^ox_si(3,j)) , ...
			    [3 1 2]);
   end
   for j = 1:size(oy_si,2)
     OMEGA(2,i,:) = OMEGA(2,i,:) + ...
                    permute( ...
                    oy_si(1,j)*(meshg.pg(1,:).^oy_si(2,j)) .* ...
                               (meshg.pg(2,:).^oy_si(3,j)) , ...
			    [3 1 2]);
   end
 end
 OMEGA = scale*OMEGA;

 for i=1:base.mk
   h(i) = figure;
   plot([0 1 0 0],[0 0 1 0],'k');
   hold on
   hh = quiver(meshg.pg(1,:),meshg.pg(2,:), ...
        permute(OMEGA(1,i,:),[2 3 1]),permute(OMEGA(2,i,:),[2 3 1]),0);
   set(hh,'color',[rand(1),rand(1),rand(1)]);
   axis equal
   title(['Basis function ',num2str(i)])
 end

 OMEGA = zeros(2,base.mk,size(pgn,2));
 for i = 1:base.mk
   ox_si = base.o_s{1,i};
   oy_si = base.o_s{2,i};
   for j = 1:size(ox_si,2)
     OMEGA(1,i,:) = OMEGA(1,i,:) + ...
                    permute( ...
                    ox_si(1,j)*(pgn(1,:).^ox_si(2,j)) .* ...
                               (pgn(2,:).^ox_si(3,j)) , ...
			    [3 1 2]);
   end
   for j = 1:size(oy_si,2)
     OMEGA(2,i,:) = OMEGA(2,i,:) + ...
                    permute( ...
                    oy_si(1,j)*(pgn(1,:).^oy_si(2,j)) .* ...
                               (pgn(2,:).^oy_si(3,j)) , ...
			    [3 1 2]);
   end
 end
 OMEGA = 0.2*scale*OMEGA;

 for i=1:base.mk
   h(base.mk+i) = figure;
   plot([0 1 0 0],[0 0 1 0],'k');
   hold on
   hh = quiver(pgn(1,:),pgn(2,:), ...
        permute(OMEGA(1,i,:),[2 3 1]),permute(OMEGA(2,i,:),[2 3 1]),0);
   set(hh,'color',[rand(1),rand(1),rand(1)]);
   axis equal
   title(['Basis function ',num2str(i)])
 end

return

